using chc 6
present CV LDAs, have testing/training if need be
graphed LDAs are all the points.
I cannot do site contrasts because of our uneven samplign scheme. The MANOVA throws out samples.
fit_full_species_man$pca_summary
## Importance of first k=7 (out of 35) components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.0736 0.7112 0.47335 0.37674 0.33980 0.30966 0.2955
## Proportion of Variance 0.4066 0.1784 0.07903 0.05006 0.04073 0.03382 0.0308
## Cumulative Proportion 0.4066 0.5850 0.66402 0.71409 0.75481 0.78864 0.8194
summary(manova(as.matrix(data[,4:cols]) ~ hostRace * sex, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## hostRace 4 0.03980 40.088 28 776.62 < 2.2e-16 ***
## sex 1 0.81535 6.956 7 215.00 1.795e-07 ***
## hostRace:sex 4 0.79391 1.833 28 776.62 0.005633 **
## Residuals 221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Reference
## Prediction Cingulata Cornivora Mendax pom Zepheria
## Cingulata 25 0 0 0 0
## Cornivora 1 5 0 0 0
## Mendax 0 0 60 5 1
## pom 0 0 1 116 0
## Zepheria 0 0 0 1 16
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 9.610390e-01 9.390180e-01 9.273307e-01 9.820323e-01 5.281385e-01
## AccuracyPValue McnemarPValue
## 1.508463e-49 NaN
Collapsed all 35 chc measurements into a value column with a corresponding chc column detailing what chc was measured. Fitted a model to evaluate the efect of host and sex on chc values. The model accounts for uneven site effects and individual effects.
All Bayesian models were created in Stan computational framework (http://mc-stan.org/) accessed with brms package (Bürkner, 2017). To improve convergence and guard against overfitting, we specified mildly informative conservative priors. We ran 4 chains with 2000 interations. We checked for convergence and found all models mixed well. (Pareto K estimates are all good!)
#1 value ~ hostRace * sex + (1|site) +(1|chc) + (1|ID)
mod = readRDS("fit_chc_model1.rds")
plot(mod)
brms::pp_check(mod, resp ="value") #great
## Using 10 posterior samples for ppc type 'dens_overlay' by default.
bayes_R2(mod, summary = TRUE)
## Estimate Est.Error Q2.5 Q97.5
## R2 0.8790622 0.0009515364 0.8771306 0.8808946
tab_model(mod)
| Â | value | |
|---|---|---|
| Predictors | Estimates | CI (95%) |
| Intercept | -0.25 | -0.49 – -0.00 |
| hostRace: Cornivora | -0.13 | -0.29 – 0.03 |
| hostRace: Mendax | 0.17 | 0.09 – 0.26 |
| hostRace: Zepheria | 0.33 | 0.23 – 0.43 |
| hostRace: pom | 0.11 | 0.04 – 0.17 |
| sex: Male | 0.03 | -0.06 – 0.11 |
| hostRaceCornivora.sexMale | -0.14 | -0.35 – 0.06 |
| hostRaceMendax.sexMale | -0.07 | -0.17 – 0.03 |
| hostRaceZepheria.sexMale | -0.01 | -0.14 – 0.12 |
| hostRacepom.sexMale | -0.14 | -0.23 – -0.05 |
| Random Effects | ||
| σ2 | 0.44 | |
| τ00 | 0.07 | |
| ICC | 0.86 | |
| N site | 8 | |
| N ID | 231 | |
| N chc | 35 | |
| Observations | 8085 | |
| Marginal R2 / Conditional R2 | 0.021 / 0.879 | |
plot_model(mod, type = "pred", terms = c("hostRace", "sex"))
plot_model(mod, bpe = "mean", bpe.style = "dot", prob.inner = 0.4, prob.outer = 0.8)
summary(manova(as.matrix(data[,4:cols]) ~ hostRace * site * sex, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## hostRace 1 0.82295 3.3194 7 108 0.003071 **
## site 1 0.85799 2.5537 7 108 0.017937 *
## sex 1 0.52503 13.9576 7 108 8.859e-13 ***
## hostRace:site 1 0.69396 6.8041 7 108 1.085e-06 ***
## hostRace:sex 1 0.95188 0.7799 7 108 0.605479
## site:sex 1 0.93727 1.0326 7 108 0.412735
## hostRace:site:sex 1 0.90238 1.6691 7 108 0.124190
## Residuals 114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Reference
## Prediction Apple Haw
## Apple 25 12
## Haw 14 71
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.786885246 0.503288443 0.703538181 0.855813669 0.680327869
## AccuracyPValue McnemarPValue
## 0.006229331 0.844519267
## Reference
## Prediction Apple_Female Apple_Male Haw_Female Haw_Male
## Apple_Female 23 2 9 1
## Apple_Male 1 1 1 2
## Haw_Female 8 0 32 3
## Haw_Male 1 3 5 30
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 7.049180e-01 5.725131e-01 6.155825e-01 7.840219e-01 3.852459e-01
## AccuracyPValue McnemarPValue
## 8.606186e-13 9.110307e-01
# no interaction because missing males at MtPleasant
summary(manova(as.matrix(data[,4:cols]) ~ sex + site, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.74708 0.67708 5 10 0.6507
## site 1 0.83749 0.38808 5 10 0.8461
## Residuals 14
## Reference
## Prediction Zepheria_Female Zepheria_Male
## Zepheria_Female 5 3
## Zepheria_Male 4 5
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.5882353 0.1793103 0.3292472 0.8155630 0.5294118
## AccuracyPValue McnemarPValue
## 0.4062810 1.0000000
## Reference
## Prediction Zepheria_Female_EastLansing
## Zepheria_Female_EastLansing 2
## Zepheria_Female_MtPleasant 0
## Zepheria_Male_EastLansing 4
## Reference
## Prediction Zepheria_Female_MtPleasant
## Zepheria_Female_EastLansing 2
## Zepheria_Female_MtPleasant 1
## Zepheria_Male_EastLansing 0
## Reference
## Prediction Zepheria_Male_EastLansing
## Zepheria_Female_EastLansing 1
## Zepheria_Female_MtPleasant 5
## Zepheria_Male_EastLansing 2
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.29411765 -0.05699482 0.10313551 0.55958272 0.47058824
## AccuracyPValue McnemarPValue
## 0.95792652 0.03207164
summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.57727 6.1024 6 50 7.541e-05 ***
## site 2 0.10716 17.1232 12 100 < 2.2e-16 ***
## sex:site 2 0.74705 1.3082 12 100 0.2258
## Residuals 55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Reference
## Prediction Mendax_Female Mendax_Male
## Mendax_Female 20 15
## Mendax_Male 13 13
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.54098361 0.07072905 0.40849889 0.66935590 0.54098361
## AccuracyPValue McnemarPValue
## 0.55241652 0.85010674
## Reference
## Prediction Mendax_Female_Fenville Mendax_Female_OtisLake
## Mendax_Female_Fenville 7 1
## Mendax_Female_OtisLake 1 5
## Mendax_Female_Sewanee 0 2
## Mendax_Male_Fenville 6 0
## Mendax_Male_OtisLake 0 1
## Mendax_Male_Sewanee 0 0
## Reference
## Prediction Mendax_Female_Sewanee Mendax_Male_Fenville
## Mendax_Female_Fenville 0 5
## Mendax_Female_OtisLake 1 0
## Mendax_Female_Sewanee 4 1
## Mendax_Male_Fenville 0 5
## Mendax_Male_OtisLake 1 1
## Mendax_Male_Sewanee 4 0
## Reference
## Prediction Mendax_Male_OtisLake Mendax_Male_Sewanee
## Mendax_Female_Fenville 0 0
## Mendax_Female_OtisLake 1 2
## Mendax_Female_Sewanee 2 3
## Mendax_Male_Fenville 0 1
## Mendax_Male_OtisLake 1 1
## Mendax_Male_Sewanee 2 3
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.409836066 0.283523654 0.285504382 0.543223627 0.229508197
## AccuracyPValue McnemarPValue
## 0.001280557 NaN
summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.78806 1.2774 4 19 0.3136059
## site 1 0.38209 7.6817 4 19 0.0007371 ***
## sex:site 1 0.70133 2.0229 4 19 0.1319094
## Residuals 22
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Reference
## Prediction Cingulata_Female Cingulata_Male
## Cingulata_Female 5 7
## Cingulata_Male 7 7
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.46153846 -0.08333333 0.26587122 0.66629178 0.53846154
## AccuracyPValue McnemarPValue
## 0.83731354 1.00000000
## Reference
## Prediction Cingulata_Female_SouthBend Cingulata_Female_Urbana
## Cingulata_Female_SouthBend 1 2
## Cingulata_Female_Urbana 2 0
## Cingulata_Male_SouthBend 4 0
## Cingulata_Male_Urbana 3 0
## Reference
## Prediction Cingulata_Male_SouthBend Cingulata_Male_Urbana
## Cingulata_Female_SouthBend 5 2
## Cingulata_Female_Urbana 4 1
## Cingulata_Male_SouthBend 0 0
## Cingulata_Male_Urbana 1 1
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.076923077 -0.243027888 0.009455391 0.251302917 0.384615385
## AccuracyPValue McnemarPValue
## 0.999943167 0.389256527
There is only one site
# no site because only sampled at one site.
summary(manova(as.matrix(data[,4:cols]) ~ sex, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.54542 0.83346 2 2 0.5454
## Residuals 3
Cant show a plot
## Reference
## Prediction Cornivora_Female Cornivora_Male
## Cornivora_Female 2 1
## Cornivora_Male 0 2
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.8000000 0.6153846 0.2835821 0.9949492 0.6000000
## AccuracyPValue McnemarPValue
## 0.3369600 1.0000000